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Abstract. In the framework of swift heavy ion - matter interaction, the thermal spike 
has proved its worth since nearly two decades. This paper deals with the necessary 
refinement of the computation due to the kind of materials involved i.e. nanomaterials 
such as multilayered systems or composite films constitued of nanocylindcrs or 
nanospheres embedded in matrix. The three dimensional computation of the thermal 
spike model is applied for the first time in layers containing spherical nanoparticles 
embedded in a silica matrix. The temperature profile calculated at each point (x,y,z) 
of the target for times up to 10~ 10 s allows a possible explanation of the particle shape 
change under irradiation with swift heavy ions having an energy of several MeV/u.m.a. 
The comparison made with the former 2D version of the code applied to cylindrical 
gold nanoparticles confirms the validity of the present 3D version. 
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1. Metallic nanoparticles synthesis with ion beam size and shape 

Since a decade, materials based on nanoparticles within a dielectric matrix have been 
elaborated with various techniques among which the use of swift heavy ions [1]. Such 
kind of materials have been found interesting because of their peculiar optical properties 
related to plasmon excitation by UV- visible radiation [2]. Potential applications are 
expected particularly in the field of optical sensors. It has been demonstrated that 
plasmon frequency depends on the particle size, on their shape (spherical, cylindrical) 
and also on the surrounding matrix [3,4]. It is then important to monitor such 
geometrical parameters with ion beam assisted method cited above. Authors use one- 
step [3, 5-7] or two-step processes [1, 8] which combine a prior ion implantation with 
the desired metallic species with post elaboration irrad iation. Since 2003 [1], several 
papers invoke the thermal spike process leading to the shaping of nanoparticles [9-11] 
and we particularly refer to the work of Ridgway et al.[12]. We have been among the first 
teams to delevop numerically this model since 1992 [13]. We have applied successfully 
this model in a wide range of materials: metals [14-16], semiconducting and insulators 
[17, 18]. In this paper, we propose the first real 3D version of the model which proves 
necessary in cases where the usual cylindrical symmetry is no longer valid: anisotropic 
thermal parameters and/or anisotropy in the geometry of the physical problem. We first 
discuss the thermal parameters that need special attention in order to give quantitative 
information instead of qualitative description. This discussion is made in relation with 
the work of Awazu et al. [19] and gives light on the choice of the thermal parameters. We 
therefore perform our simulation for the materials and ions studied in [19]: cylindrical 
gold nanoparticles embedded in silica and irradiated with Br 110 MeV. Then we focus on 
the application of the model to the nanoparticle transformation under irradiation with 
swift heavy ions and show the difference between cylindrical and spherical particles such 
as those investigated in [9-12]. Finally we provide a full description of the behaviour of 
a spherical particle due to an ion impact. 

2. The 3D thermal spike model 

The target material is considered as a two component system: the atomic lattice and the 
electrons characterized by their respective temperatures T e and T a . Starting from the set 
a coupled equations ((la) and (lb)) governing the electronic and atomic temperatures 
(T e and T a ), every physical parameter may now depend on the position f(x,y,z) in 
space. We write the energy variation dQ = CdT linked to the temperature variation 
dT of an elementary volume dV characterized by a specific heat C . dQ is due to three 
terms: i) the energy brought by the incident ion on the electrons (A e (r,t)) and on the 
atomic lattice (A a (f,t)), ii) the Fick's thermal diffusion law which defines the heat flux 
j = K ■ VT where the thermal conductivity K is a 3x3 matrix with the components 
Kij, and Hi) the energy exchange between electrons and atomic lattice proportional to 
the temperature difference g(f,t)(T e — T a ) where g(f,t) the electron phonon coupling 
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Figure 1. Geometry of the calculation box. 
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V(K e -VT e )-g(T e -T a ) + A e 
V(T a -VT a ) + g(T e -T a ) + A a 



(la) 
(lb) 



We describe the different terms in the following sections. 
2.1. Incident ion energy deposition 

The linear energy transfer from the incident ion to the target atoms (S n ) and electrons 
(S e ) were computed using SRIM code by Biersack et al.[20] and shown in the Table 
table 1. In the case of swift ions S„JS e < 1CT 2 so that we neglect the nuclear component 
A a and deal only with electronic one A e . 

We suppose that the incident ion hits the target parallel to Z-axis on (x ,y ) 
point in XF-plane (figure 1) We separate the space and time variables so that 
A e (r, t) = F(r)G(t). The function G(t) is a time pulse of few femtoseconds duration 
(to = 2 lCT 15 s), and is normalized so that f*° G(t)dt = 1. 

The function F(r) is a product F(r) = C S e (xo,yo, z) v 2 (x,y,z) W(r±). Here we 
introduce some new parameters: u(x,y,z) = pZ/A is the parameter proportional to 
the electron density (p is mass density, Z and A are atomic and mass number of target 
material); 

r±(x,y,z) = v yjjx — x ) 2 + (y — yo) 2 is the radial range (in g cm -2 ) to the ion 
track, weighted by mean electron density; 

v = J Q v [(,({}), !](-&), z] d$ where £($) = x + (x — x )$ and r]($) = y + (y — yo)$. 



Finally, W(r±) 



[l-(r ± +R 

r± (r±+Rmin) 



is known as the radial 5-ray energy 



distribution (Waligorski et al. [21]), where we substitute the range for homogeneous 
material with our weighted range r±. According to [21], R m i n , R ma x are the ranges 
corresponding to the target ionization potential E min and maximum energy of <5-electrons 
E max = 2mc 2 (3 2 /(l — (3 2 ), (3 is incident ion velocity in lightspeed units; R = k E a , 
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k = 6 l(T 6 g cm" 2 keV~ a ; a = 1.079 for /3 < 0.03 and a = 1.667 for p > 0.03; 
W(r ± > R max ) = 0. 

C is a normalizing constant, so that f F(x, y, z)dxdy = S e (z). 

The energy density transferred to the 5 electrons is plotted (in eV/atom) in two 
XY planes for two values of Z (figure 2). Hence, we see the distributions in the upper 
plane (Z = Z\) for which the ion track is outside the nanoparticle (figure 2 a)) and in 
the bottom plane [Z = Z 2 ) (figure 2 b)). Due to the different electronic densities of 
the nanoparticle and the matrix, we notice that the energy distribution is wider within 
the nanoparticle than in the matrix (e.g. see the light blue contour corresponding to a 
zone in which the energy density exceeds 10 eV/atom). The different shapes of spatial 
energy distribution proves the influence of the nanoparticle / matrix interface at this 
first step of the ion-matter interaction. 

2.2. Lattice parameters 

The lattice specific heat C a and thermal conductivity K a values are a compilation of 
data from [22-25]. However, for gold, known as a noble metal, the main contribution 
to measured thermal conductivity comes from electronic part, especially at low 
temperature. So, in order to estimate the atomic thermal conductivity, we use the 
relation K a = C a v s X/3 , where v$ is the sound speed and A is the phonon mean free 
path (of order of few lattice parameters), in the same way, C a is taken for classical 
statistical physics. The temperature dependence of C a and K a are plotted in figure 3 
and figure 4 respectively (continuous lines: blue for Au and red for SiC^). 

The latent heat of phase change, Q m , is accounted using the following algorithm 
described in the case of solid-liquid change that occurs at T = T m . For temperatures 
between T m +AT and T m +AT, we modify the temperature dependence C a {T a ) by adding 
the peak-like function C m (T a ) to C a (T a ), so that j^ m ^^ C m (T)dT = Q m , AT = 5K 
(See insert in figure 3). Since the value of C a becomes high around melting point, the 
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Figure 3. Electronic (C e in dotted 
lines) and atomic (C in straight 
lines) specific heats of Au (blue) 
and Si02 (red) as a function of 
temperature 



Figure 4. Electronic (C e in dotted 
lines) and atomic (C a in straight 
lines) thermal conductivities of Au 
(blue) and SiC>2 (red) as a function 
of temperature 



Material 


Si0 2 


Au 


Ps[g/cm 3 ] 


2.32 


19.3 


p L [g/cm 3 ] 


2.32 


17.3 


T m [K] 


1972 


1337 


T V [K] 


3223 


3130 


Qm[J/g] 


142 


63.7 


Q„[J/g] 


4715 


1645 


#[W/(cm 3 K)] 


1.25 10 13 [19] 


2.3 10 10 [15] 


S e [keV/nm] 


9.7 


29.2 


S'rtfkeV/nm] 


0.025 


0.113 



Table 1. Physical parameters of SiC>2 and Au. 



lattice temperature during melting remains practically constant. Similar modification 
is done also for vaporization account. 

2.3. Electronic parameters 

We consider the electronic system in gold as free electron gas. At low temperatures 
C e (T e ) = - k ^ ne jr; is linear with temperature according to [26]. Here T F = 

2m k B (37r 2 »-e) 2//3 is the temperature corresponding to the Fermi energy, h and are 
Planck and Boltzmann constants, m e and n e are electron mass and number density, 
respectively. The high temperature classical value C e = ^ksn e is valid for T e > -^Tp . 
C e (T e ) is presented in figure 3 (blue dotted line). 

Concerning the thermal conductivity, the relation K e = C e D e is used. The 
temperature dependence of the electronic diffusivity D e (T e ) was studied by Martynenko 
et al. [27] and used here. It was shown that the value of D e decreases from ~ 10 2 cm 2 
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s- 1 at 300K down to ~ 1 cmV 1 at T e ~ 10 4 K" 1 . So, we approximate D e (T e ) = 300^ 
at T e < T Hm , and D e = D min at T e > T Hm , with D 300 = 150 cm 2 s _1 , D min = 2 cm 2 s _1 , 
and Tu m = 300 f? 3 "" . Hence the value K e is plotted versus T e in figure 4 (blue dotted 
line). 

For the dielectrics, the free electron gas model is valid only when the temperature 
exceeds the bandgap (E g = k B T g ps 9eV leading to T g 10 5 K for Si0 2 ). We use a 
high temperature (T > T g ) specific heat value of 1 J cur 3 K) [17], which corresponds 
to hot electrons in SiC-2, whereas below this temperature, we refine the model and use 
a linear law C e (T e ). C e is plotted as function of T e in figure 4. According to the same 
reference, the high temperature thermal conductivity of SiC>2 is 2 W cm -1 K -1 . We 
refine in the same way the temperature dependence K e (T e ) using a linear law. Both 
linear laws simply account for the fact that, in a first order approximation, the number 
of electrons involved in the thermal process is proportional to the temperature T e as 
long as T e < T g . Thus, as a refinement to ref. [19], we take this into account (See 
figure 4). 

2. 4- Electron-phonon coupling constant g 

The electron-phonon coupling coupling constant g has been studied in supraconductors 
[28, 29] and the formulation of g has been successfully extended to other type of materials 
(metals [14,15,30,31], semiconductors and insulators [18,32]) . We use the values 
reported in these latter references. The value of g and K e used here are consistent with 
the ones used by [18,32] since the mean free path A (A 2 = K e /g) remains unchanged 
between the present paper and those references. 

Considering all the remarks above, we can now report in table table 1 all the 
thermodynamical parameters used for the calculations. 

3. Results and discussion 

3.1. Comparison with previous calculations: cylindrical particles 

In this section, we make the link between the present study and the work previously 
published by Awazu et al. [19]. In order to confirm the accuracy of the present 3D 
version of the thermal spike code, we have performed simulations first in bulk gold and 
silica, and then in the case of gold cylindrical nanoparticles with 10 nm and 20 nm 
diameter embedded in silica with thermodynamic parameters used in [19]. The results 
are perfectly consistent with [19] regarding the maximal temperatures reached as well 
as the heating and cooling times. We plot both electronic and atomic temperature 
evolutions as a function of time (abscissa) and radial distance to the ion path (ordinate) 
for 20 nm particle (figure 5). Like for simulation done in [19] the gold nanoparticle do 
not melt. This is due to the fact that electronic heat distributed over NP spreads out 
and is transferred to silica atoms before 10 -12 s, while this time is too short for the gold 
lattice to heat up due to the low ^-factor. 
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Figure 5. Cylindrical gold nanoparticle (diameter: 20nm). Evolution of the electronic 
and atomic temperatures T e (r, t) and T a (r, t) in Kelvin as a function of radial distance 
r from the ion path and time i.The thermodynamic parameters arc those from Awazu 
et al. [19] 



However, elongation of such nanoparticles [19], and even larger ones (40 nm [33]) 
was observed experimentally. 

We show the very different behaviour (see figure 6) if we use the Si0 2 K e {T e ) value 
(figure 4) as previously discussed (see sect, section 2.3). Up to t w 2 1CT 14 s, the 
electronic temperature for both cases is similar, but for greater times, (figure 5a), the 
temperature increase spreads out of the nanoparticle, while hot electrons are confined 
within the NP (figure 6a) so that T e remains higher than 3000K for times as long as 
3 10 -12 s inside the NP as well as in the thin SiC^/NP interface layer. Due to high g- 
factor electronic heat is efficiently transferred to the silica atoms so that the SiC>2 lattice 
gets also hot in this interface layer. Finally, due to lattice thermal diffusion, the lattice 
temperature of metal nanoparticle rises. The heating of nanoparticle atoms comes from 
the outer layers toward the core until complete NP melting(see figure 6b). This process 
takes time from 10~ 13 s to 3 10~ 12 s, then both electronic and lattice subsystems are in 
equilibrium and cool down until the NP is solidified at about 5 10 -11 s. 

3.2. Spherical gold nanoparticles 

Actually, the simulations from previous section, both ours and those described in 
[19], were performed in 2D-plane perpendicular to swift heavy ion track which avoids 
investigations as a function of depth (^-coordinate). Now, we deal with spherical 
particles embedded in a matrix so that the third coordinate is required since the 
incoming ion first penetrates the matrix and then the particle. We again compare the 
simulations resulting from the parameter set used by Awazu et al.[19] and the present 
one for two particle sizes. 
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Figure 6. Cylindrical gold nanoparticle (diameter: 20nm). Evolution of the electronic 
and atomic temperatures T e (r, t) and T a (r, t) in Kelvin as a function of radial distance 
r from the ion path and time t. The thermodynamic parameters and the thermal 
conductivities K esiC>2 and K aAu are those of the present work. 



3.2.1. 20nm gold nanoparticle: role of the thermal conductivity at the Si02 interface 
In figure 7 the maximum lattice temperature reached in 20 nm diameter spherical 
particle is shown. One can see that in case of Awazu's parameters (left hand side), 
only silica regions close to the ion path or NP surface are hot, since they are provided 
with high electronic temperature at early stage. The outside area and particle itself 
remain relatively cold. In contrary, if we account for the thermal conductivity (see 
sect, section 2.3) (right hand side), we observe very hot track in silica, and completely 
molten nanoparticle with vaporized outer layer, while outer space is almost at initial 



1337 K 

1972 K 

3130 K 
iNP edge ™=> 




Figure 7. Maximum lattice temperature reached around a gold nanoparticle (diameter 
= 20 nm). Physical parameters: a) Awazu et al. [19], b) present paper. Contour lines 
correspond to: Au melting temperature (1337K), S1O2 melting temperature (1972K), 
Au vaporizing temperature (3130K). The dotted blue line is the boundary of the 
metallic particle. 
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temperature. These results are similar to previous section, but the extra heat comes to 
polar NP regions from neigbouring silica. 

3.2.2. AOnm gold nanoparticle: influence of the gold lattice thermal conductivity 
Another peculiarity can be found for larger particles which are not completely molten 
( figure 8). Concerning the proper choice of lattice thermal conductivity of a metal, as 
we mentioned above, the measured values are not adequate since they mostly represent 
the electronic part of thermal conductivity. The molten phase on the plots is presented 
with double curves, at T m ±AT (1332K and 1342K). It can be seen that, with measured 
values of K a (figure 8a), an extended region reaches melting point (T m _ contour), though 
only a small fraction has surpassed the melting latent heat (T m+ contour). With the 
correct K a value (figure 8b), both contours coincide and the volume of complete melting 
is much larger, while the core temperature remains lower than on figure 8a. 




-40 -30 -20 -10 10 20 

distance from NP center (nm) 



-10 10 

distance from NP center (nm) 



a) 



b) 



Figure 8. Maximum lattice temperature reached in a gold particle (d — 40nm), 
using different thermal conductivities (K a ): a) K a taken from Awazu [19], b) our 
proposal (see sect, section 2.2). Contour lines correspond to the melting points: 
1332K - T m _(gold), 1342K - T rn+ (gold), 1972K - T m (Si0 2 ). 3130K is the vaporizing 
temperature of gold. 



3.3. Melting and vaporization of nanoparticles with different sizes 

The volumic molten and vaporized fraction vs. initial particle diameter is shown on 
figure 9 for both spherical and cylindrical particles with two kinds of lattice thermal 
conductivity of gold and interfacial electronic thermal conductivity. As we showed 
above, if we do not account for the interfacial region, melting is not observed even for 
20nm particles. 

The relative molten and vaporized volume is systematically higher for 3D-spherical 
particles in comparison with cylindrical ones. We can suppose that the energy deposited 
in the nanoparticle by the incident ion is then confined and evenly distributed over the 
whole volume of the particle. The value e = (S e l)/(VN a ) gives the energy density per 
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Figure 9. Molten (a) and vaporized (b) fraction vs. particle diameter 

atom received from a swift heavy ion. Here I is the ion path inside the particle, V 
and N a are particle volume and atomic concentration. For a cylinder with a radius 
R, £ cy i = while for a sphere e sp h = (in case of central ion hit). Naturally, 

molten (and vaporized) fraction is higher for spherical NP than for cylindrical one of 
the same diameter. Also, the heat exchange between metal and silica around ion spot 
plays important role. 

Complete melting is observed up to 30nm diameter. Then, if ion energy is not 
high enough to melt the whole particle, then with experimental (high) lattice thermal 
conductivity of gold the molten fraction quickly drops to zero, while using K a derived 
from the equation in section section 2.2, the molten fraction slowly decreases, remaining 
non-zero up to (70-80)nm. 

As for vaporization, there is no difference between K a Au from [19] and present work, 
since above melting point both parameters are the same. The particles with diameter 
from 10 to 15 nm are partially vaporized, while below lOnm complete vaporization is 
observed. Such small particles undergo dissolution and provide precipitated material 
for larger particles. 

3.4- Full 3D: effect of off- axis impact 

3.4-1. Temperature distribution Using the 3D capabilities to treat any problem without 
any symmetry, we focus on the effect of an ion off-axis impact (figure 10). We 
performed simulation in NP diameter range from lOnm to 80nm at different offset values: 
A = dXo/Ro = n/4, n varied from 1 to 6. 

The most interesting to see what happens when the ion hits close to the edge of 
NP. The corresponding plots are shown in figure 11 for 20nm and 40nm nanoparticles 
and A = |, 1, 1 (the central impacts for them are shown on figure 7b,figure 8b). Two 
remarks arise from the results: i) in non central impacts, the melting of the particle still 
starts from the surface, ii) the lattice temperature significantly increases only if the ion 
path intercepts the particle. 

This last point leads to the fact that, for such particles, a minimum fluence $ m i n 
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of the order of ions cm" 2 is required to transform all the particles i.e. to obtain 
a fusion at the particle surface. § m in = 8 10 10 cur 2 (respectively 3.2 10 11 cm~ 2 ) for 
particle diameter of 40nm (respectively 20nm) Refering to [12], the minimum fluences 
cited giving place to particle elongation is of the order of some 10 11 ions cm" 2 that 
means that a given particle has been hit by several ions which by statistical effect act 
as if the mean impact was centered so that the elongation eventually occurs along the 
irradiation beam direction. 

3.4-2. Energy density and molten volume figure 12 demonstrates the influence of the 
nanoparticle diameter on the deposited energy density (e) and the molten volume 
fraction (V m /Vo) for spherical and cylindrical particles, e and V m are plotted as a 
function of the ion path offset A. First, there is a difference depending on whether 
the incident ion enters NP (A < 1) or not (A > 1). In the latter case, only a small 
amount of energy is deposited due to secondary 5— electrons penetrating inside the 
particle. However this is enough to melt nanoparticles smaller than ~20nm in radius. 
Then, comparing the deposited energy density for spherical e sp h = f^Vl — A 2 and 
cylindrical particles e cy i = one can see that for A = y/E/3 ~ 3/4, e sph = e cy i- The 
e{d) curves really coincide for this latter case (figure 12a), but due to electronic heat 
leakage from the particle to the surrounding silica in the ion spot area less spherical 
particle is molten in diameter range (25-60) nm. However, starting from 60nm diameter, 
all spherical particles intercepted by ion (A < 1) demonstrate stable molten fraction 
above ~ 2%, because heated particle-Si0 2 interface keeps the temperature high enough 
to melt the metal. 

We estimate the energy density needed to melt the whole particle as e m = 
(J T ™ C a dT a + psQ m )/N a . For gold e m = 0.43 eV/at. However, comparing figure 12a 
and b, one can see that the energy necessary for complete NP melting is approximately 
twice as large as e m (from 0.8 to 1 eV/atom). This is due to heat leakage outisde the 
particle and can be taken into account for estimation with other particles and ions. The 
corresponding curves for cylindrical particles do not depend on the value of the impact 




Figure 10. Off-axis impact layout. 
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Figure 11. Maximum lattice temperature as a function of the ion path distance d 
(d = R.A) from the center of a gold nanoparticle of radius R. Left column: particle 
diameter D = 20nm; from top to bottom: A = 4,1, |. Right column: particle 
diameter D — AOnm; from top to bottom: A = 4, 1, | 
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Figure 12. Deposited energy density (a) and molten volume fraction (b) versus the 
gold nanoparticle diameter for different values of the offset A (see legend figure 11 for 
the definition of A) 



offset A provided that < A < f (i.e. when the ion strikes the particle). 



4. Conclusion 



We presented the first real 3D implementation of thermal spike model devoted to 
the simulation of swift heavy ion interaction with any kind of composite structures. 
A quantitative description of the ion shaping of metallic nanoparticles embedded in 
amorphous matrix can be done provided that the thermodynamical parameters follow 
the characteristics described hereafter. The target electrons are considered as a quasi free 
electron gas with classical specific heat and thermal diffusivity, whereas the target lattice 
is characterized with thermal coefficients based on statistical physics. The measured 
values cannot be used for metals since the thermal properties essentially depend on the 
electrons so that none of the subsystems (electrons and lattice atoms) is well described. 
For example, in metals, the electronic thermal conductivity is much higher than the 
lattice thermal conductivity. Using the 3D approach, we can make the distinction 
between cylindrical and spherical particles and show the influence of energy exchange 
on the metal/silica interface around the ion spot. Using the metallic (nanoparticle) and 
dielectric (matrix) thermal properties, it is shown that the heating of nanoparticle starts 
from the outer layer toward the core. This provides conditions for partial melting when 
there is not enough energy to melt the whole nanoparticle. Finally we demonstrate 
the influence of the impact parameter (i.e. radial distance between the ion path and 
the spherical nanoparticle center): i)a large particle (radius 60 nm) partially melts in 
the case of an ion path nearly tangent to the particle surface due to large NP/matrix 
interfacial small particle (60 nm ) can be molten even if the ion does not 

penetrate it due to heat transfer from molten silica track of about 10 nm diameter. The 
simulations made in present work essentially increase the limits of probability to melt 
metal nanoparticles irradiated with swift heavy ions and, thus, provide background for 
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understanding of ion shaping effect. 
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